Continuous Matrix Product States for Quantum Fields 
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We define matrix product states in the continuum limit, without any reference to an underlying 
lattice parameter. This allows to extend the density matrix renormalization group and variational 
matrix product state formalism to quantum field theories and continuum models in 1 spatial dimen- 
sion. We illustrate our procedure with the Lieb-Liniger model. 

PACS numbers: 



The numerical renormalization group (NRG) of Wil- 
son [l[ and the density matrix renormalization group 
(DMRG) of White revolutionized the way strongly 
correlated quantum systems can be simulated and under- 
stood. The applicability of those approaches has been 
better understood during the last 5 years by rephras- 
ing those methods in terms of matrix products states 
(MPS) d H; the success of NRG and DMRG relies on 
the fact that those MPS give a very accurate description 
of the correlations and entanglement present in ground 
states of 1-D quantum spin systems [jj @ . This insight 
led to several important extensions of DMRG, as MPS 
can also be used to describe dynamical properties Q and 
can be used as a stepping stone for constructing higher- 
dimensional analogues known as projected entangled pair 
states (PEPS) 0. 

In this paper, we show how this formalism of MPS 
can be adopted to describe quantum field theories. We 
will define a new family of states that we call continu- 
ous MPS (cMPS) that describe field theories in 1 spatial 
dimension. We will also show that cMPS can be under- 
stood as the continuous limit of standard MPS. Those 
cMPS can be used as variational states for finding ground 
states of quantum field theories, as well as to describe 
real-time dynamical features. Just as MPS capture the 
entanglement structure of low-energy states of quantum 
spin systems, the cMPS seem to capture the entangle- 
ment features of the low-energy states of quantum field 
theories. We will illustrate this on the hand of simula- 
tions that we have done on the Lieb-Liniger model 
which describes a system of bosons in a one dimension 
interacting via a delta-potential; using cMPS with a very 
low bond dimension, the ground state energy density is 
already reproduced with extremely good precision. We 
will also show how one can calculate other interesting 
physical quantities, like correlation functions or the static 
structure factor. 

Let us next define the cMPS, which is most easily done 
in the formalism of second quantization. We will con- 
sider a one-dimensional system of bosons or fermions on 
a ring of length L and associated field operators ip(x) 
with canonical commutation relations, [ip(x),ip(y)' ! }± — 
8(x — y) with < x,y < L space coordinates. A cMPS 



is defined as 



IX) = Tr a 



Ve 



f L dx[Q(x)®t+R(x)®^ (x)] 



|n) 



with Q(x),R(x) position dependent matrices of dimen- 
sion D x D that act on a Z?-dimensional auxiliary sys- 
tem, V exp the notation for the path-ordered exponen- 
tial, Tr aux the trace over the auxiliary system, and 
the vacuum state [^(x)|n) = 0]. A translational invari- 
ant state can easily be obtained by choosing Q(x) and 
R(x) independent of x, and a system with open bound- 
ary conditions can be obtained by replacing the Tr atlx 
by a left and right multiplication of the auxiliary system 
with a row and a column vector, respectively. 

As we will show later, cMPS appear very naturally as 
a continuous limit of MPS. Thus, they automatically in- 
herit all the properties of MPS, like the fact that the 
entanglement entropy of a contiguous block of bosons is 
bounded above by 21og 2 (-D). In general, the state \x) is 
a superposition of states with different particle number. 
For the case of fermions, it is easy to enforce an occu- 
pation number with a fixed parity by introducing a Z2 
symmetry by choosing Q and R block diagonal: 



Q{x) = 



Qo(x) 
Qi{x) 



R(x) = 



R (x) 
Ri(x) 



As a consequence, expectation values of the form 
{ip(xyip(y)) can be calculated without the need for in- 
troducing string-order like operators. 

To get some intuition about the structure of such 
states, it is instructive to write down explicitly 



n=Q J0<X ± 



dxi . . . dx n 4> n i>(xi) . . .-0(x„)|n), 



<...<X„<L 



where 



(fin l^oiia: 

[uq(xi, 0)Ruq(x2, xi)R . . . Ruq(L, x n )] 

and UQ(y,x) — V exp [J Q(x)dx\ . One can interpret 
uq as a free propagator, while R can be understood as 



2 



a scattering matrix that creates a physical particle. In 
general, the MPS formalism can indeed be rephrased as 
a representation of scattering events that happen in the 
physical vacuum of the interacting many-body state. 

With the help of this definition, it is straightforward 
to express the norm and expectation value of operators 
in terms of the matrices R and Q. For the sake of sim- 
plicity, we will consider a bosonic system and assume 
translational invariance. Note that for inhomogeneous 
systems one can proceed in a very similar way. Using the 
commutation relations of the field operators one readily 
finds ( X \x) = Tr (e TL ) and 



d 2 ' 



Tr [e TL (R®R)] , 



Tr 



3 T(L-x) ( R R yTx ( R R j 



dx 2 



il>{x)) = Tr[e TL ([Q,i?]® [Q,R])] , 



where T = Q®1 + 1<8)Q + R<8)R (the bar indicates 
complex conjugation). The state \x) is invariant under 
the "gauge" transformation Q -> XQX' 1 , R -> XRX- 1 
for arbitrary invertible X. This allows us to fix a gauge 
by imposing Q + + R^ R = 0, so that we can write 

Q = --R^R-iH 

where H = and R^R is diagonal. Making the trans- 
formation X®Y\a, b) — > X\a) (b\Y^ , T is transformed into 
a superoperator T (mapping matrices into matrices) , and 
we obtain that p(x) := e Tx p satisfies a master equation 
in the Lindblad form 

j-p(x) = -i[H, p(x)} + Rp(x)rf - ~ [RlR, p(x)] + . 

As a consequence, all eigenvalues of T have a non-positive 
real part, which implies that all the above quantities are 
well behaved in the thermodynamical limit L — > oo. In 
a generic case, the master equation will have a unique 
steady state p ss > 0, which can be chosen with unit 
trace. In such a case, the above expressions consider- 
able simplify in the thermodynamic limit, since (xlx) = 
Ti(p S s) = 1, 



(i>(x)H(x)) = Tr [RtR Paa ] 



{-4>{Q)^{x)^{x)^{0)) 
d 2 - 



Tr 



dx 2 



rl>(x)) 



(Re Tx (Rp ss rf)R^ 
Tr [([Q,R\y[Q,R]p ss ] 



(1) 



Other quantities can be similarly calculated. Note that 
observables defined as Fourier transforms of correlation 
functions, like the static structure factor, can be directly 
calculated in terms of the super-operator (T — ik) -1 . 

In the case of a system with open boundary conditions, 
the eigenvalues of the matrix p(x) would exactly corre- 
spond to the squares of the Schmidt coefficients when 



considering a bipartition at site x. This can in its turn 
be used to calculate the entanglement entropy of the re- 
duced density matrices defined on given intervals. Just as 
in the case of quantum spin systems, the justification for 
using cMPS should stem from the fact that an area law 
is satisfied for this entanglement entropy, eventually with 
logarithmic corrections in the case of critical systems. It 
seems indeed possible to generalize the work of Hastings 
Q (proving the area law for 1-dimensional gapped spin 
systems) to the current continuous setting [lOj. 

Let us next show how these continuous MPS can be un- 
derstood as a limit of a family of MPS. For simplicity, we 
will consider a translational invariant system of bosons 
on a ring of length L; an identical construction works for 
the fermionic case. We define a family of translational 
invariant MPS oi N = L/e modes on a discretized lat- 
tice with lattice parameter e with modes a, that obey the 

commutation relations d\,dj 



A = t + eQ 
A 1 = eR 

A n = e n R n/ n] 
dj 



Again, \Q) is the empty vacuum on which the opera- 
tors hi act (cii|fi) = 0), and we use the convention that 



(V'!) = 1- The operators ipi are defined as rescaled 
annihilation operators ; obviously, those will become the 
field operators in the limit e — > 0: 



4}] 



"'j 
e 

[i>(x),4>(y?] 



5(x - y). 



Q and R arc D x D matrices, and the scaling of the 
matrices A 1 as a function of e has been chosen such that 
this limit is well defined. The matrices A k for higher 
k have been determined by the requirement that e.g. a 
doubly occupied site yields the same physics as 2 bosons 
on 2 neighboring sites in the limit e 0. With this 
convention, the continuum limit of this MPS is equivalent 
to the continuous MPS defined before. It can also be 
checked that all divergencies in 1 /e magically cancel each 
other, such as occurring in the case of calculating the 
kinetic energy 



E, 



kin 




\x)- 
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The cancellation of the divergent terms 1/e 2 and 1/e can 
easily be proven by expanding 



$1+1- 



■0J 



[Q, R] is the only term 



as a series in e; the term [Q, R] 
independent of e. 

Let us next illustrate how these continuous MPS can 
be used as a variational ansatz for strongly correlated 
continuous theories by applying them on the Lieb-Liniger 
model Q . The Lieb-Liniger Hamiltonian describes (non- 
relativistic) bosons in 1 spatial dimension interacting via 
a contact potential: 



dx 



di[j' (x) dip(x) 
dx dx 



ctp' (x)ip' (x)ip(x)ij:(x) 



It is well known that, in the limit L — > oo, the energy den- 
sity in this system can be expressed as E/L — p 3 e(c/p) 
with p the density and e(c) the energy density of the 
system at p = 1. This scaling can also readily be un- 
derstood from the continuous MPS ansatz: for L — > oo, 
e LT remains invariant under the scaling transformation 
Q xQ and R — ► y/xR. Since density, kinetic, and in- 
teraction energy behave like R x R, [Q, R] ® [Q, R], and 
i? 2 <g)i? 2 , respectively, we have that under this transforma- 
tion p -> xp, E kin -> x 3 E kin and E lnt -> x 2 E int . Thus, 
E km {p) + 2cE int {p) = p 3 (E Hn (p = 1) + c/pE mt (p = 1)) 
giving the above scaling. The energy density e(c) can 
be determined in terms of the Bethe ansatz [9( , whereas 
other quantities whereas other quantities like correlation 
functions have been calculated using Monte Carlo meth- 



ods Ref. 11, 12 



We did a variational optimization of cMPS as a func- 
tion of the scaling parameter c (i.e., we chose p = 1). 
We carried out a simple gradient minimization of the en- 
ergy density as a function of the matrices A — iH and 
R = OD, where A is antisymmetric, O orthogonal and D 
diagonal. In Fig. 1 we have plotted e(c) for different val- 
ues of the bond dimension D, as well as the one obtained 
by Bethe ansatz; The insert shows the relative error in 
the energy as a function of D, which seems to indicate an 
exponential dependence. As a comparison, for c = 2 the 
Bethe ansatz gives e = 1.0504, the Monte Carlo method 
of Ref. 



11 



|12| gives e = 1.0518, whereas we obtain 
e = 1.1241, 1.0618, 1.0531, 1.0515, 1.0512, and 1.0508 for 
D = 2,4, ... ,12. In Fig. 2 we have determined the one- 
particle and density-density correlation functions. With 
little numerical effort we obtain results which are compa- 
rable to those of exact Monte Carlo methods. By using 
more sophisticated techniques to perform the minimiza- 
tion we believe that much more precise results can be 
obtained, and thus cMPS can be viewed as an alterna- 
tive to other existing methods 1^, 14 1 . Importantly, the 
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Figure 1: Energy density as a function of the interaction pa- 
rameter c for different values of D = 2, 4, 6, 8 (from top to 
bottom). The result for D = 8 is indistinguishable from the 
one given by the Bethe Ansatz. The insert shows the relative 
error AE = (e - e Bethe ) /e B ethe) (where E Bethe is the energy 
given by the Bethe Ansatz solution), as a function of D for 
c = 0.2,2,20 and 200 (x,*,+, and o„ respectively). We show 
the results for up to D = 10; the saturation of the accuracy 
with D = 10 is due to insufficient convergence of the results. 
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Figure 2: (a) Off-diagonal elements of the one-particle re- 
duced density operator as a function of the distance in a log- 
arithmic scale for (from top to bottom): c = 0.2, 2, , 20, 200 
and c = 2000. For reference we have also drawn a straight 
line with slope=l/2, which is the slope corresponding to the 
Tonk- Girardeau limit (c — s> oo). As it can be seen, the slope 
of the curves approaches 1/2 as c increases, (b) Two-body 
density-density correlation function for the same values of c 
(in the left, from bottom to top). For large c one can observe 
the Friedel oscillations corresponding to the Tonks- Girardeau 
limit. All the results have been calculated with D = 14. 



cMPS method does not rely on the fact that the model 
is integrable, and works equally well for non-intregrable 
models. 

Let us next comment on how to do the calculations in 
the case the translational invariance is broken. This is ob- 
viously of central importance for the simulation of atomic 
gasses in a non-homogeneous potential such as occurring 
in optical lattices; the present ansatz allows to deal with 
the full Hamiltonian as opposed to effective Hamiltoni- 
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ans such as the Bose-Hubbard model which typically ig- 
nore the potentially important effects from the higher 
Hubbard bands. In that case, one should expand the 
functionals Q{x), R{x) as a series, in such a way that a 
discrete amount of parameters characterize the state: 

p p 

Q(x) = J2 fpWQp R ( x ) = E A-W^ 



k=0 



k=Q 



Here the functions f p (x) can be chosen to correspond 
to harmonic Fourier functions in the case of a periodic 
lattice or by localized functions in the case of e.g. a har- 
monic trap. For periodic lattices, this leads to a Bloch- 
like ansatz, and it is possible to define eigenf unctions of 
the site-dependent Lindblad operator in terms of a simi- 
lar Fourier series. Similarly, it is possible to incorporate 
the MPS techniques for real-time evolution. In this case, 
the Q(x,t) and R(x,t) become both functions of space 
and time, and it is possible to write down coupled differ- 
ential equations that describe the evolution. 

Other obvious extensions include the simulation of sys- 
tems with different types of fermions and/or bosons. This 
is relevant for the case of the Hubbard type models, where 
there are 2 types of fermions per site or in the case of 
mixtures. In this case, the cMPS ansatz becomes 



Tr 

J- 1 n 



where the vj a are field operators corresponding to differ- 
ent spins (or species). Obviously, more local terms can 
be added in the exponential, such as 

S a p 0) ® (x)4>l (x) + S' a(} (x) ® ip a (x)4>l (x) . 

a/3 

Besides that, it is possible to extend this formalism to 
2-dimensional continuum systems using the formalism of 
PEPS §|. In that case, the auxiliary bond dimension has 
to be interpreted as representing an auxiliary field, and 
the judicious choice of tensors Q and R allows to develop 
a consistent formalism for describing 2+1 dimensional 
field theories fioj ]. 

In conclusion, we have introduced a new family of 
states, the cMPS, for quantum field models in 1 spatial 
dimension. They correspond to the continuum limit of 
the MPS. We have shown how one can efficiently deter- 
mine expectation values of different observables, so that 
they can be used to approximate ground state of such sys- 
tems. There are many possible extensions of the present 



work. On the one hand, one can apply the same tech- 
niques as with MPS to describe mixed states or systems 
at finite temperature, as well as higher dimensions 
On the other hand, it would be interesting to explore 
new methods for finding the matrices Q and R variation- 
ally with high bond dimension, as well as to study non- 
translationally invariant systems. Beyond that, it would 
also be interesting to substitute those matrices by opera- 
tors acting on an infinite-dimensional Hilbert space as in 
[H| in order to capture critical phenomena and to study 
relativistic quantum field theories. Finally, the cMPS 
formalism allows to construct Hamiltonians whose exact 
ground states are known, which leads to new solvable 
field theories (Toj . 

This work was supported by the EU Strep project 
QUEVADIS, the ERC grant QUERG, the FWF SFB 
grants FoQuS and ViCoM, and the DFG-Forschergruppe 
635. 



[1] 
[2] 

[3] 



[4] 



[5] 
[6] 



u. 



[8] 



[9] 
[10] 

[11] 
[12] 

[13] 
[14] 

[15] 



K.G. Wilson Rev. Mod. Phys. 47 773 (1975) 
S.R. White Phys. Rev. Lett. 69 2863 (1992) 
Schollwock Rev. Mod. Phys. 77 259 (2005) 
M. Fannes, B. Nachtergaele, R.F. Werner, Commun. 
Math. Phys. 144, 443 (1992); S. Ostlund and S. Rommer, 
Phys. Rev. Lett. 75, 3537 (1995); F. Verstraete, D. Por- 
ras, and J. I. Cirac, Phys. Rev. Lett. 93, 227205 (2004). 
F. Verstraete, J.I. Cirac, V. Murg, Adv. Phys. 57, 143 
(2008); J.I. Cirac, F. Verstraete, J. Phys. A: Math. 
Theor. 42, 504004 (2009) 

F. Verstraete and J.I. Cirac Phys. Rev. B 73 094423 
(2006) 

M.B. Hastings J. Stat. Phys. P08024 (2007) 

G. Vidal Phys. Rev. Lett. 93 040502 (2004); S.R. White 
and A.E. Feiguin Phys. Rev. Lett. 93 076401 (2004); A.J. 
Daley, C. Kollath, U. Schollwock and G. Vidal J. Stat. 
Mech.: Theor. Exp. P04005 (2004); F. Verstraete, J.J. 
Garcfa-Ripoll, and J.I. Cirac, Phys. Rev. Lett. 93 207204 
(2004) 

F. Verstraete and J.I. Cirac arXiv:cond-mat/0407066 
see also G. Sierra and M.A. Martin-Delgado, 
|cond-mat/9811170 



E. H. Lieb and W. Liniger, Phys. Rev. 130, 1605 (1963) 

F. Verstraete and J. I. Cirac, in preparation. 

G. E. Astraharchik and S. Giorgini, Phys. Rev. A 68, 
031602 (2003) 

G.E. Astraharchik and S. Giorgini, J. Phys. A 39, 1 

(2006) 

J.S. Caux. !arXiv:0908.1660l 

J. S. Caux and P. Calabrese, Phys. Rev. A 74, 031605 

(2006). 

J.I. Cirac and G. Sierra, larXiv:0911.3"029l 



